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We perform Monte Carlo simulations to study the elastic properties of the helix-coil worm-like 
chain model of alpha-helical polypeptides. In this model the secondary structure enters as a scalar 
(Ising like) variable that controls the local chain bending modulus. We compute numerically the 
bending and stretching compliances of this molecule as well as the nonlinear interaction between 
stretching and torque over a wide range of model parameters. The numerical results agree well with 
previous mean-field and perturbative calculations where they are expected to do so. The Monte 
Carlo simulations allow us to examine the response of the chain to large forces and torques where 
the perturbative approaches fail. In addition we extend our mean-field analysis by studying the 
fluctuation dominated regime at the force-induced denaturation transition. 
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I. INTRODUCTION 



The mechanical properties of semiflexible biopoly- 
mers such as F-actin are well described by the worm- 
like chain Hamiltonian 0, Q, which treats the filament 
as a one dimensional elastic continuum without inter- 
nal structure. Experiments based on single molecule 
manipulations |3, LjJ- Li Lj- LZl Lj- liJ- lii! have demonstrated 
the essential validity of this coarse-grained approach to 
the investigation of biopolymer statistics and mechanical 
properties. Understanding both scattering function of 
such semiflexible polymers in dilute solution^]] and the 
force extension relations 12] of single semiflexible poly- 
mers are among the principal successes of this coarse- 
grained description of these polymers. A simplification 
of this nature, although valid in numerous contexts, must 
breakdown under forces that are large enough to mod- 
ify the local structure of the polymer. For example the 
worm-like chain Hamiltonian ignores the double helical 
structure of DNA and the local secondary structure of 
polypeptides. This local molecular structure controls the 
bending modulus of the polymer and can be disrupted 
under applied stress. This limitation of the worm- like 
chain (WLC) model is becoming increasingly apparent as 
researchers probe the mechanical properties of biopoly- 
mers under larger applied forces where such a coupling 
between local molecular structure and chain elasticity 
plays an important role [l^. Such large forces are not 
only experimentally accessible but are also biologically 
relevant in such processes as those associated with DNA 
looping[ll[l3 and in protein conformational change. An 
example of the latter process can be found in the confor- 
mational change of the protein Calmodulin that involves 
the buckling of a single, solvent-exposed alpha-helical do- 
main upon the binding of Ca 2+ ions|16lll7|. A polymer 
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model of the alpha-helical domain that incorporates this 
nonlinear elastic response of the chain is therefore re- 
quired to explore conformation change. We expect that 
such a model (as presented here) will be broadly appli- 
cable to the investigation of protein allostery. 

Recently we0, ^3 and others 0, 113] proposed a 
minimal extension of worm-like chain called the he- 
lix/coil worm-like chain (HCWLC) to describe biopoly- 
mers with internal structure. This model couples the 
conformational degrees of freedom of the polymer back- 
bone to localized structural transitions of the constituent 
monomers by postulating that the local bending modu- 
lus of the worm-like chain depends on the local degree of 
internal structure (the helix/coil variable). For example, 
the bending modulus of an alpha-helical polypeptide de- 
pends on the local presence of secondary structure and 
the consequent hydrogen bonding that effectively stiff- 
ens the chain. The fundamental result of this coupling 
is to make both the torque and force response of the 
polymer highly nonlinear due to localized denaturation 
events (loss of local secondary structure) under mechani- 
cal stress. These denatured regions introduce more com- 
pliant elements into the chain leading to abrupt changes 
in the effective moduli of the polymer at a critical applied 
stress. 

It is important to note that the secondary structure 
variables are locally coupled along the chain since the 
alpha-helical structure is stabilized by hydrogen bonding 
between adjacent turns of the helix. These nearest neigh- 
bor interactions between the secondary structure vari- 
ables make the denaturation transition cooperative. We 
will parameterize this cooperativity by an energy scale 
in what follows referring to it as the chain cooperativity 
parameter. The worm-like chain Hamiltonian and our 
extension of it, the HCWLC are discussed in more detail 
in section ITT1 

The numerical work reported here complements our 
previous analytic calculations of the extension and bend- 
ing compliances of alpha- helical polypeptides as modelled 
within the framework of the HCWLC. For the case of ex- 
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ternally applied torques with no externally applied ten- 
sile stress, these analytic calculations have been carried 
out exactly. Both the force extension relations of the 
molecule and the nonlinear coupling of the extensional 
compliance to the state of torque, however, cannot be 
solved in closed form. A similar issue arises in the study 
of the extensional compliance of the worm- like chain |l2fl. 
Previously we examined the response of the molecule to 
tensile stress as well as to a combination of torque and 
tensile stress by two routes: (i) perturbation theory ap- 
plicable to the low force regime, and (ii) mean field the- 
ory that replaces fluctuating internal variables along the 
chain by their self-consistent, mean value. This latter 
approximation is of particular utility for the case of high 
chain cooperativity that suppresses strong fluctuations 
of the secondary structure variables. Here we provide 
a Ginzburg criterion|2l| to discuss the limits of validity 
of this approximation and use Monte Carlo simulations 
to evaluate the response of the molecule in those regimes 
where the previous mean-field or perturbative approaches 
do not apply. We also use the Monte Carlo method to 
confirm our previous calculations where such analytical 
calculations are applicable. 

This work thus explores a new range in the parameter 
space of the model, i.e. the regime of strong secondary 
structure fluctuations. Although estimates 19] show that 
at least some alpha-helical polypeptides under physiolog- 
ical conditions have high chain cooperativity and thus are 
well-approximated by the mean field description, explor- 
ing the elastic properties of a polypeptide with low chain 
cooperativity may well reflect the effects of poor solvent 
quality. In addition we will show that the model predicts 
that there is always a fluctuation-dominated regime near 
the point at which the alpha-helix denatures under ten- 
sion. 

To briefly review the phenomenology of the model, we 
find a highly nonlinear elastic response of the alpha-helix 
when subjected to bending torques and stretching forces 
as a result of the interaction between the local secondary 
structure and the conformational fluctuations. In par- 
ticular we reported a mechanical instability of the chain 
under bending akin to a "buckling instability" in which 
the segment of alpha-helical polypeptide locally melts to 
form a denatured (random coil) segment beyond a criti- 
cal angle. At this point the torque required to hold the 
chain at that particular angle drops precipitously. 

Upon application of a tensile stress to an alpha-helical 
polymer we find three regimes of response. At the lowest 
stresses, the external force pulls out the small equilib- 
rium undulations of the stiff alpha-helical polymer in a 
manner identical to the WLC. There are however, signif- 
icant differences between the force extension behaviors 
of the WLC and the HCWLC. In the WLC model the 
extensional compliance of the chain vanishes in the limit 
of large forces because each segment of the chain is in- 
extensible. The applied tensile force simply depletes the 
equilibrium population of transverse undulations along 
the chain. As these undulations vanish no additional arc 



length can be recovered by additional force; the mea- 
sured extension of the chain plateaus at its total con- 
tour length. This is not the case for the HCWLC. For 
short enough chains (to be quantified below) the stan- 
dard Marko-Siggia plateau in which additional force does 
not produce extension is replaced by what we term a 
"pseudoplateau" characterized by the fluctuations of the 
polymer into the random coil, or denatured state. Since 
the random coil sections of the chain are longer than the 
same segments in their native, alpha-helical state, the 
biasing of the fluctuations into the random coil leads to 
additional contour length of the polymer. Finally in the 
third regime the chain is completely denatured by the ap- 
plied force and reaches a true Marko-Siggia plateau upon 
further increasing the force. 

The remainder of the paper is organized as follows. 
In section [H] we introduce the alpha- helix Hamiltonian 
based on a combination of the worm-like chain and the 
helix/coil model. The details of the Monte Carlo sim- 
ulation is described in IIIII The simulation results are 
analyzed in the section IIVI We first study the radius 
of gyration of the alpha helix as a function of solvent 
quality (to be defined therein) in section HV Al We then 
examine the response of the chain to bending torques 
in section FtVB II In section FtV B 21 we study the exten- 
sional compliance of the alpha-helix before summarizing 
and discussing the results in section M 

II. THE HELIX-COIL WORM-LIKE CHAIN 
MODEL 

The worm-like chain (WLC)0, Q is the simplest 
coarse-grained model for semiflexible polymers. This 
model describes the single-chain polymer statistics in 
terms of a Hamiltonian that associates an energy cost 
with chain curvature by introducing a bending modulus 
k. In terms of a discretized chain model described by the 
set of monomeric tangent vectors ti, i = 0, . . . ,N— 1 with 
N the degree of polymerization, the WLC Hamiltonian 
may be written as 

N-l 

H WLC (ii ' U+i)] ■ (1) 

i=0 

The bending modulus k determines the thermal persis- 
tence length of the chain, the typical distance along the 
chain over which the tangent vectors decorrelate. 

The response of a single chain to extensional forces has 
been used to understand the deformational properties of 
biopolymers and their aggregates [T2I I22I |23| . To account 
for the internal degrees of freedom along the chain a new 
set of variables is needed. Workers have previously em- 
ployed the helix/coil (HC) model[2^| in order to intro- 
duce such internal state variables along the arc length of 
the chain. This model has been used to study a class of 
protein conformational transitions |25l l26| in solution and 
under tension[2^|. 
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FIG. 1: Schematic figure of an alpha-helical polypeptide and 
its representation in terms of the Ising-like secondary struc- 
ture variables (open circles for random coil segments and filled 
ones for alpha-helical ones) and the tangent vectors to the 
segments of the chain (denoted by arrows) . 



The HC model Hamiltonian, which is used to study 
these structural transitions can be reduced to its simplest 
form by assuming that the local structure of the chain can 
be described by a set of two-state variables s, = ±1, i = 
0, . . . , N. For the alpha-helical chains of current interest 
we regard these two states as the local conformation of 
the segment in its native, alpha-helical state (s = +1) 
and in a disordered, random coil state (s = —1). 

It is important to note that the elementary units of 
the chain as described by the HCWLC model are not 
the amino acid monomers but rather turns of the alpha 
helix since it is necessary to unambiguously ascribe the 
presence or absence of secondary structure to each ele- 
mentary unit of the model. In practice we expect this 
level of coarse graining to mean that each segment of the 
chain (or elementary unit) is composed of ~ 3 monomers. 

The coupling of the secondary structure variables to 
the WLC tangent vectors is effected by introducing a 
bending stiffness in the WLC Hamiltonian that depends 
on the local degree of secondary structure. We choose 



k> if s = +1 
k < if s = — 1 



(2) 



Due to the hydrogen bonding between turns of the alpha- 
helix, it is reasonable to expect that k> the bending 
modulus in the native state is significantly larger than 
k< the bending modulus of the chain in the non-native, 
disordered state. Simple estimates of this difference in 
bending moduli were computed in |19| . 

The Hamiltonian of the helix-coil worm-like chain may 
be thus written as 



N—l 



N 



H = e w /2 (1 - - h/2j2(si - 1) + 

N-l 

«(«<) [i - (u ■ , (3) 



logarithm of the chain cooperativity parameter; h rep- 
resents the free energy cost per monomer to be in the 
non-native (i.e. random-coil) state, while k> and k< are 
the bending moduli of the chain in the helix and coil 
phases as mentioned above. A pictorial representation of 
the system is shown in figure ^ 

The full Hamiltonian given by Eq.[3]has four constants 
with dimensions of energy: k>, k<, h, e w that can be fit 
from experiment. We have disregarded the twist degree 
of freedom of the molecule. Such twist degrees of free- 
dom and the coupling of twisting and stretching modes 
of these chiral molecules have been explored particu larly 
with regard to the mechanical properties of DNA[28ll29|. 
Also in the present work we evaluate the model in two di- 
mensions. Full three dimensional variants of the calcula- 
tion that incorporate torsional modes as well as the twist- 
stretch coupling are currently under investigation!^. 



III. MONTE CARLO SIMULATIONS 

We performed Monte Carlo simulations of the helix 
coil worm-like chain model using a standard Metropolis 
algorithm. Our system consists of a linear chain of N 
lattice points with two variables specified at each lattice 
position, a scalar variable Si = ±1 recording the presence 
or absence of secondary structure and an angular variable 
Qi specifying the angle the tangent vector joining two sites 
of the chain make with the mean direction. The number 
of tangent vectors in the simulation is one less than the 
total number of spins, i.e. it is iV — 1. 

The energy of the chain is specified by the Hamiltonian 
given by Eq. |3J We choose for our initial configuration a 
state having equal numbers of helix and coil segments and 
with randomized tangent vectors. There are two classes 
of Monte Carlo moves: spin flips representing attempts to 
change the local secondary structure and tangent vector 
moves allowing the chain to explore all conformations. 
To update the configuration of the chain, we flip a spin 
at random Sj — > — and calculate the energy difference 
between the final (v) and the initial configurations (/i). 
The acceptance probability of this move is given by 



A(ji -> v) 



I 



otherwise 



E„ > 



(4) 



where e w is the free energy cost of a domain wall in the 
secondary structure sequence and is the negative natural 



where [i and v represent the initial and final state of the 
system respectively. 

In the next sweep of the lattice we pick Qi at random, 
attempt a Monte Carlo update of an angular variable 
Qi — ► Qi ± AQ and accept or reject it according to the 
Metropolis algorithm outlined above. We find that choos- 
ing AQ = 0.01 allows us to equilibrate the chain over 
reasonable timescales while being small enough so that 
the discretization of chain conformations does not signif- 
icantly affect the numerical results. This latter point was 
checked by a comparison to runs with even smaller val- 
ues of the angular variable updates. We study systems of 
sizes ranging from N — 10 to N — 100. We see a small 
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but finite system-size dependence in the sharpness of the 
transitions of the underlying spin variables. 

Finally, because of their discrete nature the spin vari- 
ables equilibrate much faster than the angular variables. 
In all the runs we checked that we had performed a suffi- 
cient number of Monte Carlo steps in order to equilibrate 
the largest correlated structures in the system that have 
the slowest dynamics, i.e. the equilibration time r eq was 
always taken to be greater than 4r cor - the largest corre- 
lation time in the system. Due to the wide separation of 
the time scales for the equilibration of the helix-coil and 
angular variables, a more efficient Monte Carlo scheme 
could be developed by updating the angular variables 
more often than the secondary structure variables when 
deep in the ordered phase of these secondary structure 
variables. We did not pursue such improvements of the 
efficiency of the code there. 

We benchmarked our code by comparing our results 
to the known equilibrium properties of the Ising model 
and the Kratky-Porod model in the limit where we have 
artificially frozen either the tangent vectors or the spin 
variables respectively. In these cases the HCWLC model 
reduces to these well-studied cases. 

As an example of these checks of the numerics we first 
examine the WLC limit of the model by fixing the sec- 
ondary structure variables Sj = +1 for all i. In this 
limit, the HCWLC Hamiltonian is given by Eq.QJup to a 
trivial constant. The correlation function for this WLC 
Hamiltonian can be calculated exactly 2] . For our tan- 
gent vectors in two dimensions we find that 

(U ■ ij) ~ exp[-|j - i\/l p ], (5) 

where l p is the persistence length of the chain that is 
given in terms of the bending modulus by 

l P = -log[—^], (6) 

This length is the distance along the chain over which 
tangent vectors decorrelate. We have equilibrated the 
chain over 2 x 10 4 time steps (corresponding to about 
500 lattice sweeps) before making our measurements. We 
show in panel a) of the figure Fig|2] the dependence of 
the persistence length upon the bending modulus. The 
dashed line shows the agreement to the prediction of 
Eq.O In panel b) of the same figure we examine the Ising 
limit of the model in which we fix the tangent vectors to 
lie along the i-axis. Here we plot the zero-temperature 
phase transition of the helix-coil (Ising) variables as a 
function of the solvent quality h (magnetic field in the 
Ising language). As predicted for the one-dimensional 
Ising model, the system jumps spontaneously to the he- 
lix phase upon making h an arbitrarily small positive 
quantity. 
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FIG. 2: a) Variation of the persistence length as a function 
of bending modulus of the WLC with k measured in units of 
ksT . b) The Ising limit of the model showing the magnetiza- 
tion curve: M vs h at ksT — 10 -4 . Both the figures are for 
an N = 20 chain. The error bars are the size of the symbols. 

IV. RESULTS 

A. The radius of gyration 

We now turn to the full HCWLC model and discuss 
first the radius of gyration[Tlj. To do so we note that 
the projection of the polymer arc length along the aver- 
age tangent vector of a segment, i.e. the effective length 
of the segment depends on the state of secondary struc- 
ture. To account for this aspect of the coarse-grained 
HCWLC polymer model we define a segment length that 
is a function of the secondary structure variable s n via 

where, as the notation suggests, 7< < 7>. The length 
of a segment increases when it loses its alpha-helical sec- 
ondary structure. We return to a discussion of reasonable 
numerical estimates of these values in the conclusions. 

The separation vector between the i th and j th segments 
along the chain is given by 

j-i 

Rij = ^7( s n)4, (8) 

n—i 

where j(s n ) is the length of the rfi 1 segment measured 
along its mean chain tangent as defined above. 

To compute the radius of gyration in our simulations 
we first evaluate the center of mass of the polypeptide. 
This is given by the expression 

1 N j 

j=i »=i 

The radius of gyration is then evaluated by computing 
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FIG. 3: Radius of gyration of the HCWLC as a function of the 
free energy cost per segment to transform to the random coil, 
nonnative state: h obtained from Monte Carlo simulations 
(plotted using black dots) as compared with the theory [l9j|. 
In this curve «> = 100, ft< = 1, N = 10, and e w = 10. All 
the energy scales are measured in fcsT and the error bars are 
smaller than the symbols in the plot. 



the average 

r g = - i) )■ (10) 

j=i \i=i / 

We plot in Fig. 01 the radius of gyration as a function 
of h for e w = 10, and bending moduli k> = 100, and 
k< = 1 for an N — 10 sized chain obtained using Monte 
Carlo simulations (black dots). These data are in excel- 
lent agreement with our analytic calculations (solid line). 
The energy scale h is the free energy cost of a segment 
being in its nonnative (i.e. random coil) state and so 
variations of h mimic changing the solvent quality as is 
done in the denaturation of proteins by the addition of 
e.g. Guanidinium[3]| . A second effect of solvent quality 
is the change in the effective excluded volume interaction 
between segments of the polymer. This is not included 
here but, since the chain never approaches a random coil 
configuration in these simulations we expect that the ab- 
sence of self-avoidance plays a minimal role in the above 
results. 

Similar behavior is seen in the dependence of the 
squared averaged end-to-end vector of the chain. For 
high chain cooperativities and h chosen so that the chain 
is in either the all helix or all coil phase we find this 
quantity has the following expected scaling form for a 
worm-like chain: AR 2 = N 2 j 2 g(-^-), where we have in- 
troduced the function: g(x) = 2(exp[— x] — l + x)/x 2 and 
7 is either 7> or 7< for the all helix or all coil chains re- 
spectively. The behavior of the HCWLC is similar to the 
WLC but with a persistence length that interpolates be- 
tween that of the helix and the random coil. This results 
also agrees with previous analytic calculations 19J. 



B. Mechanical Properties 

1. Torque response 

To begin our numerical exploration of the mechanical 
properties of the model, we consider the response of the 
chain to externally applied torques in thermal equilib- 
rium. These torques act to constrain the tangent vectors 
of the ends of the molecule while not applying tensile 
stress. Since the restricted partition function of the chain 
with arbitrarily fixed end tangents can be computed ex- 
actly, we expect the results of the simulations to agree 
with the previous work for all values of the model pa- 
rameters. The response of the chain to tensile stress or 
to the combination of tensile stress and bending torques, 
however, cannot be computed analytically in closed form; 
it is there that the complementary utility of the Monte 
Carlo approach will be seen. 

To study the torque response of the chain, we hold the 
first tangent vector fixed along the i-axis and constrain 
the last chain tangent to make a fixed angle ip with re- 
spect to the same axis. We numerically evaluate the con- 
straining torque r(ip) by computing the derivative of the 
free energy with respect to the angle ip. In our simula- 
tions this is effected by directly calculating the thermal 
average: 

r(tp) = (/s(sjv-i) sin[^ - 6 N -i]). (11) 

The results, which corroborate the analytic calcula- 
tions are plotted in figure^] At small values of the bend- 
ing angle ip, there is a linear dependence of the constrain- 
ing torque on ip. The alpha-helix bends like a flexible, 
elastic rod. At a certain critical angle ip*, however, the 
constraining torque reaches a maximum and then drops 
precipitously for angles ip > ip* as shown in figure 01 
This dramatic collapse of the chain's rigidity is akin to 
the buckling instability of a macroscopic tube such as a 
drinking straw. The mode of the localized failure though 
is quite different. In the case of current interest, the 
failure is caused by the localized disruption of the sec- 
ondary structure. The breaking of the hydrogen bonds 
at this denatured site introduces a weak link allowing 
the molecule to bend at a lower torque. At ip = ip*, M, 
the fraction of the chain in the nonnative state, abruptly 
jumps to 0(1/ N) demonstrating that, within the model, 
the buckling failure is due to the creation of a single 
random coil segment along the chain that provides a re- 
gion of greatly reduced bending stiffness. The size of the 
created random-coil section will remain on the order of 
Nk < /k > so for a large difference in bending moduli be- 
tween the native and nonnative states of the chain, these 
"weak links" generically occupy a small fraction of the 
polymer. For instance, in the example shown in figure 01 
there is only one weak link created. 

The above behavior is seen in the parameter range (e w , 
h) consistent with the chain being in an all helix state 
in thermal equilibrium in the absence of applied torque. 
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FIG. 4: Numerical torque vs. angle data (filled points) for 
e w ~ 5.2, h = 8, K> = 100, and «< = 1 compared with 
previous exact calculations [13,14] (solid line) for an N = 10 
chain. The torque is measured in units of ksT. The error 
bars are of the same size as the symbols in the plot. 



If the equilibrium system is in a mixed helix-coil phase 
such a behavior is not observed. Instead we see a mono- 
tonic increase in the bending torque for all end angles 
less than 7r as expected for an elastic rod. At ip — tt 
the torque measured in thermal equilibrium is identically 
zero since one is averaging a signed quantity and the sta- 
tistical weight of positive and negative torques becomes 
equal. The average of the squared torque, however, van- 
ishes only at tp = modulo 2ir. Since it is reasonable to 
suppose that the chain bends in the plane defined by first 
and last constrained chain tangents, one expects that our 
two-dimensional results accurately captures the buckling 
instability of the chain. We expect that the full three 
dimensional calculation would generate at least qualita- 
tively similar results. 

We now turn to the numerical study of the extensional 
compliance of the chain. These Monte Carlo simulations 
complete the analysis of extensional compliance of the 
model by accessing parameter regimes where neither of 
our previous approximations are valid. In particular we 
look at the extensional compliance where the mean field 
approximation fails and we look at the modification of 
the extensional compliance due to applied torque. The 
existence of this torque/stretch coupling is due to the 
inherent nonlinearity of the model. 



2. Extensional compliance 

In the presence of a tensile force F, the Hamiltonian 
of the HCWLC may be written as 



N 



(12) 



i=Q 



where Ho is the HCWLC Hamiltonian in the absence of 
externally applied forces as shown in Eq.|3J In our Monte 



Carlo simulations we apply a force of magnitude F and 
then equilibrate the chain as described above. We then 
calculate the quantity 



JV 



(13) 



i=l 



This result gives the mean length of the chain as a func- 
tion of the externally applied force. Similar measure- 
ments are made for the case where the initial and final 
chain tangents are constrained: Q\ = and 9n = ift- 
These latter simulations allow us to explore the nonlinear 
dependence of the extensional compliance upon applied 
torque. In such simulations we compute 



JV 



L{F^)= 



(14) 



In all the extension measurements the force is always 
collinear with the initial chain tangent, t\ = x. 

We plot in figure \5\ the force-extension curve of the 
HCWLC model with two sets of parameter values. In 
both force extension measurements we do not constrain 
the final chain tangents and numerically determine the 
function defined in Eq. In the first case (upper panel 
of the figure) we have selected the HCWLC model param- 
eters such that the polymer is in a fluctuation-dominated 
regime where the mean-field approximation is not valid. 
In the lower panel we have tuned those parameters to 
suppress secondary structure fluctuations. 

In upper panel we take e w = 0.5 and h = 1.0 to en- 
hance the secondary structure fluctuations. The bending 
moduli are chosen to be k> = 4 and k< = 2 respectively. 
Here we do not expect the mean field theory (shown as 
a solid line in figure [SJ) to fit the simulation data. The 
mean field theory indeed fails to describe the force ex- 
tension curve adequately. In particular we note that the 
most significant discrepancy between the mean-field the- 
ory and the data occurs at low forces. In the limit of 
high forces, the tension in the polymer acts as an order- 
ing field suppressing the fluctuations that the mean-field 
theory ignores. Moreover, we note that the intermedi- 
ate plateau ( "pseudoplateau" ) in the graph of extension 
vs. applied force vanishes for this parameter regime. We 
return to the vanishing of the pseudoplateau below. In 
fact in the limit of low chain cooperativity, the polymer is 
typically found in a mixed helix/coil phase and we obtain 
force extension curves similar to those the WLC results 
of Marko and Siggiafl^. 

In the limit of high chain cooperativity (simulation 
data and mean field prediction shown in the lower panel 
of figure the data agree well with our previous mean- 
field calculations using the HCWLC model as expected. 
In this case we have taken e w = 8 to enforce high chain 
cooperativity and suppress secondary structure fluctua- 
tions. We have taken the remaining parameters to be 
h = 1.5, k> = 100 and k < = I. 
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FIG. 5: Force vs. extension Monte Carlo data (points) shown 
in the low (upper panel) and high (lower panel) cooperativity 
limit. In each the mean length of the chain normalized by the 
maximum chain length A r 7> is plotted as a function of the 
applied force normalized by the length of a helix segment 7> . 
For the upper panel the parameters are e w = 0.5, h = 1.0, 
«> = 4.0, k< = 2.0, and N = 20. The lower panel the 
parameters are e w = 8, and h = 1.5, and k> = 100, k< = 1, 
and N = 10. In both graphs the mean field calculation [l9| 
is shown as a solid line. 



Returning to the appearance of the pseudoplateau in 
the upper panel of figure we recount the explanation 
of its existence. For small applied forces, a molecule in 
the all-helical phase having large stiffness approaches its 

maximum length iV7< as 1/ \ 



^2^. The WLC predicts 
the vanishing of the extcnsional compliance as the chain 
reaches its maximal length; this would produce a plateau 
at these intermediate forces. This plateau, however, is 
not flat in the HCWLC due to the fact that an increase 
in the applied force enhances the fluctuations into the 
longer, random coil phase of the segments, which we re- 
fer to as the pseudoplateau. The end of this plateau is 
marked by a sharp lengthening transition at a value of 
force F + ss e \~t h where A7 = 7> — 7< . Here we observe 



A7 



the force- induced denaturation of helical domains. 

Upon closer examination one notes that the appear- 
ance of this intermediate-force pseudoplateau is actually 
a finite-size effect in the model. The width of this pseu- 
doplateau is controlled by two forces, a force at onset F_ 
and a maximum force F + at which point the abrupt de- 
naturation transition occurs. These forces are marked on 
figure The lower force F_ is determined by the begin- 
ning of the first Marko-Siggia plateau where the ampli- 
tude of transverse thermal undulations of the all alpha- 
helical polymer have been significantly reduced. Because 
that reduction has an algebraic dependence on applied 
force, it is not possible to unambiguously select a critical 
force marking the onset of the plateau. It appears rea- 
sonable, however, to insist that plateau has been reached 
when dL/dF (1/t<) <C 1, i.e. the when the incremen- 
tal extension of the chain measured in monomer lengths 
5L /7< is small for a change in force SF ~ l/7< set by the 



thermal energy (k^T — 1) divided by the same monomer 
size. In this case we find the force associated with the 
plateau onset F- ~ N 2 ^ 3 k > 1 ^ 3 j^ 1 grows with the length 
of the chain. Longer chains require larger forces to suffi- 
ciently pull out the equilibrium population of transverse 
undulations and thereby reach the plateau regime. 

The high force end of the plateau is determined by bal- 
ancing the free energy cost associated with transforming 
a segment from alpha helix to random coil with the work 
done by the external force during that transformation: 
h + e w ~ A7F.4. . On the left hand side is the free energy 
cost of making one domain wall and converting one seg- 
ment to the thermodynamically unfavored random coil 
state. The right hand side of the same relation gives the 
work done by the external force as one segment length- 
ens by A7. The existence of the pseudoplateau requires 
that F_ < F + ; the applied force must straighten out the 
alpha helix at tensions small enough such that the force- 
induced helix-to-coil transition does not occur. This final 
inequality demands that 



N < 



fh 



\7< A 7 



3/2 



(15) 



Typical alpha helical polypeptides are quite short TV ~ 
O(10) and stiff k> ~ O(10 2 ) so that one expects to ob- 
serve this pseudoplateau behavior quite generally. We see 
that for the parameters used to create the upper panel 
of figure [3] the criterion for the presence of the pseudo- 
plateau is not met and in our Monte Carlo simulations we 
indeed observe no pseudoplateau for these values. When 
the criterion is met as in the case shown in the lower panel 
of figure [5] the mean field example, the pseudoplateau is 
evident. 

We now consider the force extension measurements 
made while constraining the initial and final chain tan- 
gents. In this manner we compute the mean extension as 
a function of applied force as well as the imposed total 
angular bend on the polymer. This quantity is defined 
by Eq. [21 In the regime where the persistence length of 
the chain in both the helix and the coil states is smaller 
than the chain itself, these constraining torques have es- 
sentially no impact on the extensional compliance of the 
chain. The effect of the angular constraints is limited 
to boundary regions within a persistence length of the 
constrained ends. For stiffer chains having a persistence 
length longer than the chain itself, there is a significant 
nonlinear coupling between these constraining torques 
and the extensional compliance of the polymer at low 
to moderate forces. 

Such effects are shown in figure where we plot two 
force extension curves for a HCWLC having parameters 
k> = 40, k< = 2, e w = 8 and h = 1.5. In both cases 
the initial chain tangent is constrained to lie along the 
x-axis, the direction of the tensile force. In one case 
(triangles) the final chain tangent is collinear with the 
initial chain tangent, while in the second case (circles) 
the final chain tangent makes an angle ip — ir/2 with 
the initial one. These simulations numerically evaluate 
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the thermal average of the chain's extension as defined 
in Ea. HI with ip = 0,tt/2. 

At low applied force we see that, due to the nonlin- 
ear coupling of applied torque and force, the chain has 
a greater extensional compliance in the bent state. The 
-0 = 7r/2 chain reaches its pseudoplateau more rapidly 
with increasing force. Once the chain denatures under 
the tensile stress causing the persistence length to be- 
come shorter than the chain itself, the coupling between 
the extensional compliance and end tangent constraint 
disappears. As one approaches arbitrarily high forces 
the effect of the constrained end tangent itself becomes 
localized at the last segment of chain regardless of the 
bending moduli of the polymer so that its effect becomes 
negligible. In the limiting case of infinite force there re- 
mains only the effect of the constrained end itself, which 
decreases the length of the chain by 7>(1 — cosi/j). In 
the strongly fluctuating secondary structure regime, we 
have not been able to extract any systematic information 
about the nonlinear coupling of force and torque for low 
forces. 

The numerical calculation presented here does not di- 
rectly probe the force/torque coupling since we have mea- 
sured the extensional compliance in an ensemble of chains 
with fixed end tangents. We have directly measured 
the function defined by Eq. O Clearly, the constraint 
torques applied to the end have a complicated depen- 
dence on the extensional stress. In order to quantita- 
tively probe the force-torque coupling it is necessary to 
work in a fixed torque ensemble rather than a fixed angle 
one. The results presented here may, however, directly 
apply to protein mechanics where allosteric rotations of 
the final tangents of an alpha helical domain can directly 
modify the extensional elasticity of the helix thus affect- 
ing both its mean length and fluctuation spectrum. 

Returning to the observed breakdown of the mean field 
theory for the force extension curves, we ask can one de- 
velop a simple criterion to predict the region of HCWLC 
parameter space where our analytic mean field analysis 
will hold? Indeed from the energy fluctuations observed 
in the Monte Carlo simulation one can directly map the 
boundaries of the fluctuation-dominated regime. It is 
clearly desirable to develop a criterion that can be writ- 
ten directly in terms of the model parameters so that one 
may know immediately where the simpler analytic theory 
applies. 

We observe that if the polymer were straight the effect 
of the force on the remaining secondary structure vari- 
ables would be to simply shift the effective free energy 
cost per segment to destroy the secondary structure from 
h to h e s — h — FAj; the cost of denaturing the segment 
is h but the net work done by the chain upon extension 
under the external force F is — FA7. We recall that the 
validity of the mean field theory requires only the sup- 
pression of secondary structure fluctuations. We then ex- 
pect that mean field theory to fail where h e g = h — FA^f 
is small, which for reasonable values of h requires signif- 
icant forces. Assuming that these forces have quenched 
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FIG. 6: Force vs. extension curves for HCWLC chain with 
constrained end tangents. The force vs. extension is shown 
for chains where the first tangent lies along the direction of 
the pulling force while the second tangent makes an angle 
of if) = (triangles) or ip = tt/2 (circles) with respect to 
the initial chain tangent. The parameters of the chain are: 
k> = 40 k < = 2, e„ = 10, h = 1, and N = 20. Note 
the extensional compliance depends on the state of the final 
chain tangent. The inset shows that the chain extension in the 
direction of applied force saturates at high forces to different 
values as a result of the constrained end tangent. 



most of the contour fluctuations of the chain, we take 
as an approximate criterion for the breakdown of mean 
field theory the Ginzburg criterion for the one-dimension 
helix-coil model with an effective field h c g. Thus mean 
field theory is expected to hold when 



< sf > 



< Si >' 



< 1. 



(16) 



Neglecting boundary effects to restore the translational 
invariance along the chain and taking TV large enough 
so that we may consider only the Ai, larger of the two 
eigenvalues of the transfer matrix, this condition can be 
written as 



Ai(9 2 Ai/9/i 2 ) 



< 1 



(17) 



(dXx/dh) 2 

where the larger eigenvalue, computed in is given by 
Ai = (l + e-'^) + 

- e ~' lcff ) 2 + 4exp[-2e w - h cB ], (18) 

where /i e ff is defined above. 

We determined the boundary marking the limits of va- 
lidity of the mean field theory in the parameter space 
spanned by e w and h c s is shown in Figd To the right 
of the dots mean field theory holds while in the region 
to the left of this division the system is dominated by 
large secondary structure fluctuations. Two traces of to- 
tal energy as a function of Monte Carlo time represen- 
tative of each region are shown as insets. In figure [7| 
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FIG. 7: A map of the e w — h c g parameter space of the 
helix coil degrees of freedom showing the region of validity 
of the mean field approximation (MF) as determined by the 
Ginzburg criterion. The insets show the energy fluctuations 
over Monte Carlo steps in numerical simulations of the two 
regimes. 



the filled points represent the boundary at which the en- 
ergy fluctuations reach one part in 10 3 of the mean. The 
Ginzburg criterion defined by Eq. El gi ye the dashed 
line when the left hand side of Eq. El is a ls° taken to 
be 10~ 3 . The correspondence of the Ginzburg criterion 
and the Monte Carlo fluctuation data demonstrates the 
utility of this measure of the validity of the mean field ap- 
proximation. The Ginzburg criterion as described above 
correctly distinguishes these two regimes. From that fig- 
ure it is clear that mean field theory holds in the limit of 
high chain cooperativity as expected. We note, however, 
that since the force-induced denaturation of the chain re- 
quires h c g — > we expect that at this transition there 
will be significant fluctuation effects. Such effects have 
not been fully explored. 



V. CONCLUSIONS 

We have numerically explored the nonlinear elastic re- 
sponse of the helix/coil worm- like chain under applied 
forces and torques. Our Monte Carlo simulations of the 
extensional compliance of the polymer agree satisfacto- 
rily with previous analytic calculations in two limits: (i) 
for small forces our results are consistent with pertur- 
bative calculations, and (ii) where mean field theory is 
expected to hold as determined by a Ginzburg criterion, 
we find that our mean field analysis is consistent with 
the Monte Carlo data. It should also be noted that the 
response of the polymer to applied torque agrees with the 
analytic calculations; this is to be expected as the force- 
free partition function of the system can be computed 
exactly in closed form. 

The appropriate energy scales h and e w to describe typ- 
ical alpha-helical polypeptides under physiological condi- 
tions can be estimated to be h ss 1.5 and e w ss 7 [32, 331. 



For such parameters and polypeptides of no more than 
N ~ O(10) turns it appears from our Ginzburg crite- 
rion that the chain is typically well-described by the pre- 
vious mean field approximation. Nevertheless in order 
to understand the statistical mechanics of force-induced 
denaturation of long polypeptides where A7 F ~ h or 
the denaturation via changes in solvent quality so that 
F = and h, e w — > + one must study the fluctuation 
dominated regime that is accessible via Monte Carlo. In 
this fluctuation-dominated regime we find that the ex- 
tensional behavior is well described by the Marko Sig- 
gia description of the WLC with a persistence length 
that interpolates between its value in the helix and coil 
states. In that regime the pseudo-plateau is absent. Such 
Monte Carlo calculations can, in effect, probe the un- 
folding transition state for such force-induced unfolding 
within the current model. By further developing the 
model and studying these transition states one may be 
able to assess the role of chain cooperativity in the un- 
folding process. A better understanding of this point 
should help to elucidate constant velocity, single-protein 
unfolding experiments. [3 H3- 

Moreover the Monte Carlo calculations provide a non- 
perturbative approach to investigating the nonlinear cou- 
pling between response of the chain to combinations of 
applied tensile forces and bending torques. Previously we 
had studied such nonlinear couplings only in the limit of 
small forces. Our new results discussed above allow us 
to explore the high force regime and reveal that at high 
forces torque plays little role in determining the exten- 
sional compliance of the molecule. However at low forces 
and for stiff chains such that the persistence length of 
the helical segments is comparable to, or larger than the 
length of the chain, the torque plays an important role. 
Such high force couplings may play a role in understand- 
ing protein conformational change. 

Taken in combination with the previous analytical 
treatments, these Monte Carlo calculations provide a 
nearly complete description of the mechanical proper- 
ties of the HCWLC model for alpha-helical polypeptides. 
It remains only to develop a full, three-dimensional de- 
scription of the polymer that incorporates the torsional 
degrees of freedom to the problem as well as to address 
the role of chemical heterogeneity in order to develop a 
complete and rather general equilibrium description of 
these biopolymers. 
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